﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: GaussianDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;
using BOL.Maths.Functions;
using BOL.Linq.Descriptive;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents a Gaussian distribution with mean, mu and standard deviation, sigma.
    /// <see>Mathematica: ref/NormalDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/Normal_distribution</see>
    /// </summary>
    public sealed class GaussianDistribution : IUnivariateDistribution<double>, IEquatable<GaussianDistribution>
    {
        #region Private variables

        private double _mu; // mean
        private double _sigma; // std deviation

        #endregion

        #region Public properties

        public static GaussianDistribution StandardDistribution { get { return new GaussianDistribution(0, 1); } }

        /// <summary>Gets or sets the mean of a Gaussian distribution.</summary>
        public double Mu { get { return _mu; } set { _mu = value; } }

        /// <summary>Gets or sets the std deviation of a Gaussian distribution.</summary>
        public double Sigma
        {
            get { return _sigma; }
            set
            {
                if (value <= 0)
                    throw new ArgumentOutOfRangeException("value", "sigma is expected to be a positive real number.");
                _sigma = value;
            }
        }

        /// <summary>Gets the range of a Gaussian distribution.</summary>
        public IRange<double> Domain { get { return new Range<double>(double.NegativeInfinity, double.PositiveInfinity); } }

        /// <summary>Gets the mean of a Gaussian distribution.</summary>
        public double Mean { get { return _mu; } }

        /// <summary>Gets the median of a Gaussian distribution.</summary>
        public double Median { get { return _mu; } }

        /// <summary>Gets the mode of a Gaussian distribution.</summary>
        public double Mode { get { return _mu; } }

        /// <summary>Gets the variance of a Gaussian distribution.</summary>
        public double Variance { get { return Math.Pow(_sigma, 2); } }

        /// <summary>Gets the precision parameter of a Gaussian distribution.</summary>
        public double Precision { get { return 1 / Variance; } }

        /// <summary>Gets the variance of a Gaussian distribution.</summary>
        public double Skewness { get { return 0; } }

        /// <summary>Gets the kurtosis of a Gaussian distribution.</summary>
        public double Kurtosis { get { return 0; } }

        /// <summary>Gets the entropy of a geometric distribution.</summary>
        public double Entropy { get { return Math.Log(2 * Math.PI * Math.Exp(1) * _sigma * _sigma); } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates a Gaussian distribution.
        /// </summary>
        /// <param name="mu">mean</param>
        /// <param name="sigma">stdandard deviation</param>
        public GaussianDistribution(double mu, double sigma)
        {
            if (sigma <= 0)
                throw new ArgumentOutOfRangeException("sigma", "sigma is expected to be a positive real number.");

            _mu = mu;
            _sigma = sigma;
        }

        /// <summary>
        /// The distribution with mu = 0 and sigma = 1 is called standard normal.
        /// </summary>
        public GaussianDistribution() : this(0, 1) { }

        #endregion

        #region ICloneable implementation

        public GaussianDistribution Clone()
        {
            return new GaussianDistribution(_mu, _sigma);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region Public methods

        /// <summary>
        /// Returns the Gaussian density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>Gaussian density</returns>
        public double Pdf(double value)
        {
            return Math.Exp(-0.5 * Math.Pow(value - _mu, 2) / Math.Pow(_sigma, 2)) / (Math.Sqrt(2 * Math.PI) * _sigma);
        }

        /// <summary>
        /// Returns the cummulative Gaussian density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative Gaussian density</returns>
        public double Cdf(double value)
        {
            return 0.5 * (DistributionFunctions.Erf(value / Math.Sqrt(2)) + 1);
        }

        /// <summary>
        /// Returns a random value from the distribution.
        /// <see>http://people.csail.mit.edu/milch/blog/index.html</see>
        /// <see>The implementation uses the Box-_muller transformation [G.E.P. Box and M.E. _muller (1958) "A note on the generation of random normal deviates".  Ann. Math. Stat 29:610-611].  See also http://www.cs.princeton.edu/introcs/26function/MyMath.java.html</see>
        /// </summary>
        /// <returns>random value</returns>
        public double Sample(Random r)
        {
            return _mu + (_sigma * Math.Sin(2 * Math.PI * r.NextDouble()) * Math.Sqrt((-2 * Math.Log(r.NextDouble()))));
        }

        /// <summary>
        /// Returns the value of the standard normal distribution for given probability.
        /// Need to implement non-static method for generic normal distribution.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public double Quantile(double p)
        {
            const double c0 = 2.515517;
            const double c1 = 0.802853;
            const double c2 = 0.010328;
            const double d1 = 1.432788;
            const double d2 = 0.189269;
            const double d3 = 0.001308;
            var sign = -1.0;

            if (p > .5)
            {
                sign = 1.0;
                p = 1 - p;
            }

            var t = Math.Sqrt(-2 * Math.Log(p));

            return sign * (t - (c0 + t * (c1 + t * c2)) / (1.0 + t * (d1 + t * (d2 + t * d3))));

            //// Coefficients in rational approximations
            //double[] a = new double[6] { -3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00 };
            //double[] b = new double[5] { -5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01 };
            //double[] c = new double[6] { -7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00 };
            //double[] d = new double[4] { 7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00 };

            //// Define break-points.
            //double plow = 0.02425;
            //double phigh = 1 - plow;

            //double q;
            //// Rational approximation for lower region:
            //if (p < plow)
            //{
            //    q = Math.Sqrt(-2 * Math.Log(p));
            //    return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
            //                                    ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);
            //}

            //// Rational approximation for upper region:
            //if (phigh < p)
            //{
            //    q = Math.Sqrt(-2 * Math.Log(1 - p));
            //    return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
            //                                           ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);
            //}

            //// Rational approximation for central region:
            //q = p - 0.5;
            //double r = q * q;
            //return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q /
            //                         (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1);
        }

        /// <summary>
        /// Estimates parameters of the chi square distribution using maximum likelihood.
        /// <see>http://en.wikipedia.org/wiki/Normal_distribution#Estimation_of_parameters</see>
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<double> source)
        {
            if (source == null)
                throw new ArgumentNullException("source");

            var list = source.ToList();
            Mu = list.Average();
            Sigma = list.StandardDeviation();
        }

        /// <summary>
        /// Estimates parameters of the chi square distribution using Bayesian methods.
        /// </summary>
        public void BayesianEstimate(IEnumerable<double> source, params dynamic[] priorParameters)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<GaussianDistribution> implementation

        public bool Equals(GaussianDistribution other)
        {
            return _mu.Equals(other._mu) & _sigma.Equals(other._sigma);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return Mu.GetHashCode() ^ Sigma.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is GaussianDistribution))
                throw new InvalidCastException("The 'other' argument is not a GaussianDistribution object.");
            
            return Equals(other as GaussianDistribution);
        }

        public override string ToString()
        {
            return String.Format("Normal(mu = {0:G4}, sigma = {1:G4})", _mu, _sigma);
        }

        #endregion
    }
}
